3.1.1 Daten einlesen

Die Strukturdaten sind vom Bundeswahlleiter zu beziehen; über diesen Link kommt man zu den Daten (CSV-Format).

Die Variablennamen sind hier erklärt.

d_str_file <- "https://bundeswahlleiter.de/dam/jcr/b1d3fc4f-17eb-455f-a01c-a0bf32135c5d/btw21_strukturdaten.csv"

d_str <- read_delim(d_str_file,
                    delim = ";", 
                    escape_double = FALSE,
                    locale = locale(decimal_mark = ",",
                                    grouping_mark = "."),
                    trim_ws = TRUE,
                    skip = 8) 

Hier sind die Namen der Spalten:

names(d_str)

##  [1] "Land"                                                                                                      
##  [2] "Wahlkreis-Nr."                                                                                             
##  [3] "Wahlkreis-Name"                                                                                            
##  [4] "Gemeinden am 31.12.2019 (Anzahl)"                                                                          
##  [5] "Fläche am 31.12.2019 (km²)"                                                                                
##  [6] "Bevölkerung am 31.12.2019 - Insgesamt (in 1000)"                                                           
##  [7] "Bevölkerung am 31.12.2019 - Deutsche (in 1000)"                                                            
##  [8] "Bevölkerung am 31.12.2019 - Ausländer/-innen (%)"                                                          
##  [9] "Bevölkerungsdichte am 31.12.2019 (EW je km²)"                                                              
## [10] "Zu- (+) bzw. Abnahme (-) der Bevölkerung 2019 - Geburtensaldo (je 1000 EW)"                                
## [11] "Zu- (+) bzw. Abnahme (-) der Bevölkerung 2019 - Wanderungssaldo (je 1000 EW)"                              
## [12] "Alter von ... bis ... Jahren am 31.12.2019 - unter 18 (%)"                                                 
## [13] "Alter von ... bis ... Jahren am 31.12.2019 - 18-24 (%)"                                                    
## [14] "Alter von ... bis ... Jahren am 31.12.2019 - 25-34 (%)"                                                    
## [15] "Alter von ... bis ... Jahren am 31.12.2019 - 35-59 (%)"                                                    
## [16] "Alter von ... bis ... Jahren am 31.12.2019 - 60-74 (%)"                                                    
## [17] "Alter von ... bis ... Jahren am 31.12.2019 - 75 und mehr (%)"                                              
## [18] "Bodenfläche nach Art der tatsächlichen Nutzung am 31.12.2019 - Siedlung und Verkehr (%)"                   
## [19] "Bodenfläche nach Art der tatsächlichen Nutzung am 31.12.2019 - Vegetation und Gewässer (%)"                
## [20] "Fertiggestellte Wohnungen 2019 (je 1000 EW)"                                                               
## [21] "Bestand an Wohnungen am 31.12.2019 - insgesamt (je 1000 EW)"                                               
## [22] "Wohnfläche am 31.12.2019 (je Wohnung)"                                                                     
## [23] "Wohnfläche am 31.12.2019 (je EW)"                                                                          
## [24] "PKW-Bestand am 01.01.2020 - PKW insgesamt (je 1000 EW)"                                                    
## [25] "PKW-Bestand am 01.01.2020 - PKW mit Elektro- oder Hybrid-Antrieb (%)"                                      
## [26] "Unternehmensregister 2018 - Unternehmen insgesamt (je 1000 EW)"                                            
## [27] "Unternehmensregister 2018 - Handwerksunternehmen (je 1000 EW)"                                             
## [28] "Schulabgänger/-innen beruflicher Schulen 2019"                                                             
## [29] "Schulabgänger/-innen allgemeinbildender Schulen 2019 - insgesamt ohne Externe (je 1000 EW)"                
## [30] "Schulabgänger/-innen allgemeinbildender Schulen 2019 - ohne Hauptschulabschluss (%)"                       
## [31] "Schulabgänger/-innen allgemeinbildender Schulen 2019 - mit Hauptschulabschluss (%)"                        
## [32] "Schulabgänger/-innen allgemeinbildender Schulen 2019 - mit mittlerem Schulabschluss (%)"                   
## [33] "Schulabgänger/-innen allgemeinblldender Schulen 2019 - mit allgemeiner und Fachhochschulreife (%)"         
## [34] "Kindertagesbetreuung am 01.03.2020 - Betreute Kinder unter 3 Jahre (Betreuungsquote)"                      
## [35] "Kindertagesbetreuung am 01.03.2020 - Betreute Kinder 3 bis unter 6 Jahre (Betreuungsquote)"                
## [36] "Verfügbares Einkommen der privaten Haushalte 2018 (EUR je EW)"                                             
## [37] "Bruttoinlandsprodukt 2018 (EUR je EW)"                                                                     
## [38] "Sozialversicherungspflichtig Beschäftigte am 30.06.2020 - insgesamt (je 1000 EW)"                          
## [39] "Sozialversicherungspflichtig Beschäftigte am 30.06.2020 - Land- und Forstwirtschaft, Fischerei (%)"        
## [40] "Sozialversicherungspflichtig Beschäftigte am 30.06.2020 - Produzierendes Gewerbe (%)"                      
## [41] "Sozialversicherungspflichtig Beschäftigte am 30.06.2020 - Handel, Gastgewerbe, Verkehr (%)"                
## [42] "Sozialversicherungspflichtig Beschäftigte am 30.06.2020 - Öffentliche und private Dienstleister (%)"       
## [43] "Sozialversicherungspflichtig Beschäftigte am 30.06.2020 - Übrige Dienstleister und \"\"ohne Angabe\"\" (%)"
## [44] "Empfänger/-innen von Leistungen nach SGB II  Oktober 2020 -  insgesamt (je 1000 EW)"                       
## [45] "Empfänger/-innen von Leistungen nach SGB II  Oktober 2020 -  nicht erwerbsfähige Hilfebedürftige (%)"      
## [46] "Empfänger/-innen von Leistungen nach SGB II  Oktober 2020 -  Ausländer/-innen (%)"                         
## [47] "Arbeitslosenquote Februar 2021 - insgesamt"                                                                
## [48] "Arbeitslosenquote Februar 2021 - Männer"                                                                   
## [49] "Arbeitslosenquote Februar 2021 - Frauen"                                                                   
## [50] "Arbeitslosenquote Februar 2021 - 15 bis 24 Jahre"                                                          
## [51] "Arbeitslosenquote Februar 2021 - 55 bis 64 Jahre"                                                          
## [52] "Fußnoten"

Es ist vielleicht praktisch, die Spaltennamen für spätere Verwendung in einer Textdatei abzuspeichern:

d_str_names <-
  tibble(
    var_name = names(d_str),
    ) %>% 
  mutate(id = row_number())

write_csv(d_str_names, "objects/d_str_names.csv")

3.1.2 Daten aufbereiten

Die Spaltennamen sind etwas unhandlich. Formulieren wir lieber prägnanter:

names(d_str) <- paste0("V",1:ncol(d_str))

d_str2 <-
  d_str %>% 
  select(state = V1,
         area_nr = V2,
         area_name = V3,
         for_prop = V8,
         pop_density = V9,
         pop_move = V11,
         income = V36,
         unemp = V47) 

Sichern wir diese Daten in eine Datei:

write_csv(d_str2, file = "objects/d_str2.csv")

Oder als Excel-Exitstrategie:

writexl::write_xlsx(d_str2, path = "objects/d_str2.xlsx")

3.2.1 Daten einlesen

Die Daten sind vom Bundeswahlleiter zu beziehen. Unter diesem Link kommt man direkt zur CSV-Datei.

Eine Erklärung zu den Variablen findet sich hier.

elec_results_file <- "https://www.bundeswahlleiter.de/bundestagswahlen/2021/ergebnisse/opendata/csv/kerg2.csv"

elec_results <- read_delim(elec_results_file,
                    delim = ";", 
                    escape_double = FALSE,
                    locale = locale(decimal_mark = ",",
                                    grouping_mark = "."),
                    trim_ws = TRUE,
                    skip = 9
                    ) 

Alternativ kann man die Datei selber herunterladen und im Verzeichnis der Rmd-Datei abspeichern. Man kann als CSV-Datei oder als XLSX-Datei abspeichern und dann in R importieren. Etwas nervig ist, dass das deutsche Excel (im Standard) keine Standard-CSV erstellt.

Hier nur bruchstückhaft dargestellt; so können Sie die Datendatei importieren:

elec_results <- read_csv("kerg2_aufbereitet.csv")
elec_results <- read_csv2("kerg2_aufbereitet.csv")  # Deutsches Excel!
elec_results <- readxl::read_xlsx("kerg2_aufbereitet.xslx")

kerg2_aufbereitet.csv meint, Sie haben die ersten paar Zeilen selber gelöscht und insgesamt dafür gesorgt, dass es eine schöne, maschinenfreundliche Datendatei ist.

Eine CSV-Datei, die mit dem deutschen Excel erstellt wurde, können Sie mit read_csv2(dateiname) einlesen.

Achtung: Excel schreibt manchmal in Ihre CSV-Datei auch ohne dass Sie auf Speichern klicken! Fragen Sie mich nicht, warum Excel sich diese Freiheit nimmt.

3.2.2 Daten aufbereiten

Konzentrieren wir uns auf die Zweitstimme, da die bei der BTW die entscheidende ist.

elec_results2 <- 
  elec_results %>% 
  select(Gebietsart, Gebietsnummer, Gebietsname, UegGebietsart, UegGebietsnummer, Gruppenart, Gruppenname, Stimme, Prozent, DiffProzentPkt) %>% 
  filter(Gruppenname == "AfD") %>% 
  filter(Stimme == 2)

Im anderen Datensatz wird “Bundesgebiet” mit “Deutschland” bezeichnet. Um die beiden Datensätze im Weiteren zusammenfügen zu können, sollten jeweils die gleiche Bezeichnung (für das gleiche “Ding”, sprich Gesamtdeutschland) verwenden werden:

elec_results2 <- 
  elec_results2 %>% 
  mutate(Gebietsname = ifelse(Gebietsname == "Bundesgebiet",
                              "Deutschland",
                              Gebietsname))

Sichern wir diese Daten in eine Datei:

write_csv(elec_results2, file = "objects/elec_results2.csv")

Der Befehl könnte voraussetzen, dass der Ordner objects existiert. Alternativ können Sie auch einfach ins aktuelle Verzeichnis schreiben, das ist einfacher (aber nicht so aufgeräumt):

write_csv(elec_results2, file = "elec_results2.csv")

4.1.1 … über alle Bundesländer

Achtung! Alle Bundesländer werden in der folgenden Analyse gleich gewichtet!

elec_results2 %>% 
  filter(Gebietsart == "Land") %>% 
  summarise(AfD_mean = mean(Prozent))

## # A tibble: 1 × 1
##   AfD_mean
##      <dbl>
## 1     12.1

4.1.2 … über Deutschland

elec_results2 %>% 
  filter(Gebietsart == "Bund") %>% 
  summarise(AfD_mean = mean(Prozent)) 

## # A tibble: 1 × 1
##   AfD_mean
##      <dbl>
## 1     10.3

4.1.3 … über Wahlkreise

elec_results2 %>% 
  filter(Gebietsart == "Wahlkreis") %>% 
  summarise(AfD_mean = mean(Prozent, na.rm = T))

## # A tibble: 1 × 1
##   AfD_mean
##      <dbl>
## 1     10.5

7.1.1 Ausländeranteil

d2 %>% 
  filter(Gebietsart == "Wahlkreis") %>% 
  select(Prozent, for_prop) %>% 
  ggplot() +
  aes(x = for_prop, y = Prozent) +
  geom_point(alpha = .7) +
  geom_smooth()

Interessant! Ein einfacher linearer Trend liegt nicht vor. Vielleicht sehen wir eher zwei, unterschiedliche Cluster?

Cluster 1 ist geprägt von hohem Ausländeranteil und Cluster 2 von geringem. In beiden Clustern ist der Zusammenhang negativ. Allerdings ist dieser Zusammenhang deutlich stärker ausgeprägt für Cluster 1.

7.1.2 Ausländeranteil nach Bundesland

Wenn wir diese Analyse aufteilen nach Bundesländern, sehen wir vielleicht klarer.

d2 %>% 
  filter(Gebietsart == "Wahlkreis") %>% 
  select(Prozent, for_prop, state) %>% 
  ggplot() +
  aes(x = for_prop, y = Prozent, color = state) +
  geom_point(alpha = .7) +
  geom_smooth(method = "lm") +
  facet_wrap(~ state) +
  scale_y_continuous(limits = c(0, 40)) +
  theme(legend.position = "none")

Es scheint sich in jedem Bundesland ein negativer, linearer Zusammenhang zu zeigen. Es sei denn, der AfD-Anteil ist sehr gering, dann wird der Zusammenhang schwach, also eine Bodeneffekt.

8.4.1 Modell berechnen

tic()
m1 <- stan_glm(Prozent_z ~ unemp_z,
               data = d3,
               refresh = 0)
toc()

## 1.978 sec elapsed

8.4.2 Ergebnis

m1

## stan_glm
##  family:       gaussian [identity]
##  formula:      Prozent_z ~ unemp_z
##  observations: 315
##  predictors:   2
## ------
##             Median MAD_SD
## (Intercept) 0.0    0.1   
## unemp_z     0.0    0.1   
## 
## Auxiliary parameter(s):
##       Median MAD_SD
## sigma 1.0    0.0   
## 
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanreg

Oder kürzer:

coef(m1)

##   (Intercept)       unemp_z 
## -8.329159e-05  3.733159e-02

Wie man sieht, spielt die Arbeitslosigkeit keine Rolle zur Erklärung des AfD-Wahlerfolgs.

plot(m1)

Man kann sich auch andere Darstellungen der Posteriori-Verteilung zeigen lassen:

mcmc_areas_ridges(m1,
                  pars = "unemp_z")

Offensichtlich kann es nicht ausgeschlossen werden, dass die Arbeitslosigkeit unabhängig vom AfD-Wahlerfolg ist.

m1_r2 <- bayes_R2(m1)

8.4.3 Modellprüfung

pp_check(m1)

tibble(
  resid = resid(m1),
  pred = predict(m1)
) %>% 
  ggplot() +
  aes(x = pred, y = resid) +
  geom_hline(yintercept = 0, linetype = "dotted") +
  geom_point() +
  geom_smooth()

m1_loo_r2 <- loo_R2(m1)
median(m1_loo_r2)

## [1] -0.003957106

8.4.4 Prioris

prior_summary(m1)

## Priors for model 'm1' 
## ------
## Intercept (after predictors centered)
##   Specified prior:
##     ~ normal(location = 0.00014, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 0.00014, scale = 2.5)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = 0, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 0, scale = 2.5)
## 
## Auxiliary (sigma)
##   Specified prior:
##     ~ exponential(rate = 1)
##   Adjusted prior:
##     ~ exponential(rate = 1)
## ------
## See help('prior_summary.stanreg') for more details

tic()
m1_prior_pv <- stan_glm(Prozent_z ~ unemp_z,
                        data = d3,
                        prior_PD = TRUE,
                        refresh = 0)
toc()

## 1.619 sec elapsed

Hinweise zu prior_PD finden sich hier.

coef(m1_prior_pv)

## (Intercept)     unemp_z 
##  0.04892458  0.02111651

plot(m1_prior_pv)

Die Prioris sehen vernünftig aus.

Vergleich der Priori-Verteilung mit der Post-Verteilung:

posterior_vs_prior(m1)

8.4.5 Fazit

m1_results <-
  tibble(
    model_id = "m1",
    priors = "standard",
    preds = "unemp_z",
    r2 = median(m1_r2),
    r2_loo = median(m1_loo_r2),
    coefs = list(coef(m1))
  )

results <-
  results %>% 
  bind_rows(m1_results)

8.5.1 Modell berechnen (m2)

tic()
m2 <- 
  stan_glm(Prozent_z ~ for_prop_z,
           data = d3,
           refresh = 0
           )
toc()

## 1.986 sec elapsed

8.5.2 Ergebnis

m2

## stan_glm
##  family:       gaussian [identity]
##  formula:      Prozent_z ~ for_prop_z
##  observations: 315
##  predictors:   2
## ------
##             Median MAD_SD
## (Intercept)  0.0    0.0  
## for_prop_z  -0.6    0.0  
## 
## Auxiliary parameter(s):
##       Median MAD_SD
## sigma 0.8    0.0   
## 
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanreg

Modellparameter visualisieren:

plot(m2)

mcmc_areas(m2,
           pars = "for_prop_z")

Der Ausländeranteil spielt laut unserem Modell eine große Rolle.

m2_r2 <- bayes_R2(m2)
median(m2_r2)

## [1] 0.333135

8.5.3 Modellprüfung

pp_check(m2)

m2_loo_r2 <- loo_R2(m2)
median(m2_loo_r2)

## [1] 0.3276684

8.5.4 Prioris

tic()
m2_prior_pd <- 
  stan_glm(Prozent_z ~ for_prop_z,
           data = d3,
           prior_PD = TRUE,
           refresh = 0
  )
toc()

## 1.305 sec elapsed

m2_prior_pd

## stan_glm
##  family:       gaussian [identity]
##  formula:      Prozent_z ~ for_prop_z
##  observations: 315
##  predictors:   2
## ------
##             Median MAD_SD
## (Intercept) 0.0    2.4   
## for_prop_z  0.0    2.5   
## 
## Auxiliary parameter(s):
##       Median MAD_SD
## sigma 0.7    0.7   
## 
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanreg

plot(m2_prior_pd,
      pars = "for_prop_z")

8.5.5 Fazit

m2_results <-
  tibble(
    model_id = "m2",
    priors = "standard",
    preds = "for_prop_z",
    r2 = median(m2_r2),
    r2_loo = median(m2_loo_r2),
    coefs = list(coef(m2))
  )

results <-
  results %>% 
  bind_rows(m2_results)

8.6.1 Modell berechnen (m3)

tic()
m3 <- 
  stan_glm(Prozent_z ~ for_prop_z + unemp_z,
           data = d3,
           refresh = 0
           )
toc()

## 2.062 sec elapsed

8.6.2 Ergebnis

m3

## stan_glm
##  family:       gaussian [identity]
##  formula:      Prozent_z ~ for_prop_z + unemp_z
##  observations: 315
##  predictors:   3
## ------
##             Median MAD_SD
## (Intercept)  0.0    0.0  
## for_prop_z  -0.6    0.0  
## unemp_z      0.2    0.0  
## 
## Auxiliary parameter(s):
##       Median MAD_SD
## sigma 0.8    0.0   
## 
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanreg

Interessant. Zusätzlich zum Ausländeranteil spielt die Arbeitslosigkeit jetzt doch eine Rolle: Es findet sich jetzt eine positive Assoziation mit der AV.

plot(m3, plotfun = "mcmc_dens")

Details zu plot() im Zusammenhang von rstanarm finden sich hier.

m3_r2 <- bayes_R2(m3)
median(m3_r2)

## [1] 0.3643387

m3_loo_r2 <- loo_R2(m3)
median(m3_loo_r2)

## [1] 0.3544175

8.6.3 Modellprüfung

pp_check(m3)

Hm, auch dieses Modell erklärt die Daten nicht gut.

8.6.4 Fazit

m3_results <-
  tibble(
    model_id = "m3",
    priors = "standard",
    preds = "unemp_z + for_prop_z",
    r2 = median(m2_r2),
    r2_loo = median(m2_loo_r2),
    coefs = list(coef(m2))
  )

results <-
  results %>% 
  bind_rows(m3_results)

8.7.1 Modell berechnen (m4)

tic()
m4 <- 
  stan_glm(Prozent_z ~ state,
           data = d3,
           refresh = 0
           )
toc()

## 2.198 sec elapsed

8.7.2 Ergebnis

m4

## stan_glm
##  family:       gaussian [identity]
##  formula:      Prozent_z ~ state
##  observations: 315
##  predictors:   16
## ------
##                             Median MAD_SD
## (Intercept)                 -0.2    0.1  
## stateBayern                 -0.1    0.1  
## stateBerlin                 -0.2    0.1  
## stateBrandenburg             1.5    0.2  
## stateBremen                 -0.4    0.3  
## stateHamburg                -0.8    0.2  
## stateHessen                 -0.1    0.1  
## stateMecklenburg-Vorpommern  1.4    0.2  
## stateNiedersachsen          -0.4    0.1  
## stateNordrhein-Westfalen    -0.4    0.1  
## stateRheinland-Pfalz        -0.1    0.1  
## stateSaarland                0.1    0.2  
## stateSachsen                 2.6    0.1  
## stateSachsen-Anhalt          1.7    0.2  
## stateSchleswig-Holstein     -0.5    0.2  
## stateThüringen               2.4    0.2  
## 
## Auxiliary parameter(s):
##       Median MAD_SD
## sigma 0.5    0.0   
## 
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanreg

plot(m4, regex_pars = "^state")

Sachsen und Thüringen haben die stärksten (absoluten und positiven) Koeffizienten.

m4_r2 <- bayes_R2(m4)
m4_r2_loo <- loo_R2(m4)

8.8.1 Modell berechnen (m5)

tic()
m5 <- 
  stan_glm(Prozent_z ~ Sachsen_Thueringen,
           data = d4,
           refresh = 0
           )
toc()

## 1.383 sec elapsed

8.8.2 Ergebnisse

m5

## stan_glm
##  family:       gaussian [identity]
##  formula:      Prozent_z ~ Sachsen_Thueringen
##  observations: 315
##  predictors:   2
## ------
##                    Median MAD_SD
## (Intercept)        -0.2    0.0  
## Sachsen_Thueringen  2.6    0.1  
## 
## Auxiliary parameter(s):
##       Median MAD_SD
## sigma 0.7    0.0   
## 
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanreg

plot(m5, pars = "Sachsen_Thueringen")

Das ist ein sehr starker Effekt.

m5_r2 <- bayes_R2(m5)
m5_r2_loo <- loo_R2(m5)

8.9.1 Modell berechnen (m6)

tic()
m6 <- 
  stan_glm(Prozent_z ~ state + for_prop_z + unemp_z,
           data = d3,
           refresh = 0
           )
toc()

## 2.074 sec elapsed

8.9.2 Ergebnisse

m6

## stan_glm
##  family:       gaussian [identity]
##  formula:      Prozent_z ~ state + for_prop_z + unemp_z
##  observations: 315
##  predictors:   18
## ------
##                             Median MAD_SD
## (Intercept)                  0.3    0.1  
## stateBayern                 -0.2    0.1  
## stateBerlin                 -0.6    0.2  
## stateBrandenburg             0.6    0.2  
## stateBremen                 -1.0    0.3  
## stateHamburg                -1.1    0.2  
## stateHessen                 -0.2    0.1  
## stateMecklenburg-Vorpommern  0.3    0.2  
## stateNiedersachsen          -0.9    0.1  
## stateNordrhein-Westfalen    -0.9    0.1  
## stateRheinland-Pfalz        -0.4    0.1  
## stateSaarland               -0.5    0.2  
## stateSachsen                 1.7    0.2  
## stateSachsen-Anhalt          0.7    0.2  
## stateSchleswig-Holstein     -1.1    0.2  
## stateThüringen               1.6    0.2  
## for_prop_z                  -0.3    0.0  
## unemp_z                      0.2    0.0  
## 
## Auxiliary parameter(s):
##       Median MAD_SD
## sigma 0.4    0.0   
## 
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanreg

m6_r2 <- bayes_R2(m6)
m6_r2_loo <- loo_R2(m6)
median(m6_r2_loo)

## [1] 0.8220622

8.9.3 Modellprüfung

pp_check(m6)

Das Modell erklärt die Daten gut.

8.10.1 Modell berechnen (m7)

tic()
m7 <- 
  stan_glm(Prozent_z ~ Sachsen_Thueringen + for_prop_z + unemp_z,
           data = d4,
           refresh = 0
           )
toc()

## 1.547 sec elapsed

8.10.2 Ergebnisse

m7

## stan_glm
##  family:       gaussian [identity]
##  formula:      Prozent_z ~ Sachsen_Thueringen + for_prop_z + unemp_z
##  observations: 315
##  predictors:   4
## ------
##                    Median MAD_SD
## (Intercept)        -0.2    0.0  
## Sachsen_Thueringen  2.1    0.1  
## for_prop_z         -0.4    0.0  
## unemp_z             0.1    0.0  
## 
## Auxiliary parameter(s):
##       Median MAD_SD
## sigma 0.6    0.0   
## 
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanreg

m7_r2 <- bayes_R2(m7)
m7_r2_loo <- loo_R2(m7)
median(m7_r2_loo)

## [1] 0.627569

8.11.1 Modell berechnen (m8)

tic()
m8 <- 
  stan_glm(Prozent_z ~ Sachsen_Thueringen,
           data = d4,
           refresh = 0
           )
toc()

## 1.479 sec elapsed

8.11.2 Ergebnisse

m8

## stan_glm
##  family:       gaussian [identity]
##  formula:      Prozent_z ~ Sachsen_Thueringen
##  observations: 315
##  predictors:   2
## ------
##                    Median MAD_SD
## (Intercept)        -0.2    0.0  
## Sachsen_Thueringen  2.6    0.1  
## 
## Auxiliary parameter(s):
##       Median MAD_SD
## sigma 0.7    0.0   
## 
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanreg

m8_r2 <- bayes_R2(m8)
m8_r2_loo <- loo_R2(m8)
median(m8_r2_loo)

## [1] 0.5002692

8.12.1 Modell berechnen (m8)

tic()
m9 <- 
  stan_glm(Prozent_z ~ Sachsen_Thueringen + for_prop_z,
           data = d4,
           refresh = 0
  )
toc()

## 1.41 sec elapsed

8.12.2 Ergebnisse

m9

## stan_glm
##  family:       gaussian [identity]
##  formula:      Prozent_z ~ Sachsen_Thueringen + for_prop_z
##  observations: 315
##  predictors:   3
## ------
##                    Median MAD_SD
## (Intercept)        -0.2    0.0  
## Sachsen_Thueringen  2.1    0.1  
## for_prop_z         -0.4    0.0  
## 
## Auxiliary parameter(s):
##       Median MAD_SD
## sigma 0.6    0.0   
## 
## ------
## * For help interpreting the printed output see ?print.stanreg
## * For info on the priors used see ?prior_summary.stanreg

m9_r2 <- bayes_R2(m9)
m9_r2_loo <- loo_R2(m9)
median(m9_r2_loo)

## [1] 0.616266

8.12.3 Visualisierung

Es lohnt sich vielleicht, dieses Modell zu visualisieren.

m9_draws <-
  m9 %>% 
  as_tibble() %>% 
  rename(intercept = `(Intercept)`) %>% 
  dplyr::select(-sigma)

ggplot(d4) + 
  aes(x = for_prop_z, y = Prozent_z) +
  # Ungewissheit der Regressionsgeraden:
  geom_abline(
    aes(intercept = intercept, slope = for_prop_z),
    data = sample_n(m9_draws, 100),
    color = "grey60",
    alpha = .15
  ) +
  # Mediane Schätzwerte für "die" Regressionsgerade:
  geom_abline(
    intercept = median(m9_draws$intercept),
    slope = median(m9_draws$for_prop_z),
    size = 1,
    color = "blue"
  ) +
  geom_point()